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Abstract. Weak gravitational lensing is considered to be one of the most powerful tools to study the mass and the 
mass distribution of galaxy clusters. However, the mass-sheet degeneracy transformation has limited its success. 
We present a novel method for a cluster mass reconstruction which combines weak and strong lensing information 
on common scales and can, as a consequence, break the mass-sheet degeneracy. We extend the weak lensing 
formalism to the inner parts of the cluster and combine it with the constraints from multiple image systems. We 
demonstrate the feasibility of the method with simulations, finding an excellent agreement between the input and 
reconstructed mass also on scales within and beyond the Einstein radius. Using a single multiple image system 
and photometric redshift information of the background sources used for weak and strong lensing analysis, we find 
that we are effectively able to break the mass-sheet degeneracy, therefore removing one of the main limitations 
on cluster mass estimates. We conclude that with high resolution (e.g. HST) imaging data the method can more 
accurately reconstruct cluster masses and their profiles than currently existing lensing techniques. 
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1. Introduction 

Clusters of galaxies have long been recognised as excellent 
laboratories for many cosmological applications. An espe- 
cially important diagnostic is their number density as a 
function of mass and redshift. This can only be measured 
if reliable mass estimates of the observed clusters can be 
obtained. In addition, in the framework of the ACDM cos- 
mological model, the dark mat ter distribution in cl usters 
likely follows the NFW profile l)Navarro et al.lll997|b 

Weak gravitational lensing is one of the most pow- 
erful tools currently available for studying the mass dis- 
tribution of clusters of galaxies. The first weak lens- 
ing de tection in clusters has been made by iTvson et alJ 
How ever, it was o nly after the pioneering work 
bv lKaiser fe Sauiresl l)l993[) that the field began to flour- 
ish, and since then many clu ster mass reconstructions 
have been carried out (see e.g. [Clowe fc Schneiderll200ll 
l2002t iGavazzi et aTll2004 lLombardi et alJl2005|) . A dis- 
agreement occurring in some cases between the cluster 
mass estimated from the weak/strong lensing measure- 
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ments with X-rays is still not well understood, although 
several s cenarios have been proposed to resolv e this issue 
(see e.g. lAllenllilfflsllEttori fc Lombardill2003h . 

In the absence of the redshift information, the main 
limitation for a precise we ak lensing mass estimat e is 
the mass-sheet degeneracy Schneider fc Seitdll99.il . If 
the redshifts of background sources and/or lens are not 
known, the transformation of the surface mass-density 
K — > k' = Art + (1 — A), where A ^ is an arbitrary 
constant, leaves the expec tation value o f measu red image 
ellipticities unchanged. In iBradac et alJ l)2004bl) we show 
that this degeneracy can be lifted using information on in- 
dividual source redshifts, however only if the weak lensing 
reconstruction is extended to the critical parts of the clus- 
ter. Strong lensing is affected by this transformation as 
well. Namely, the mass-sheet degeneracy does not change 
the image positions (since the source position is not an 
observable) and flux ratios and therefore can not be bro- 
ken if a single redshift multiple-image system is used. The 
mass-sheet degeneracy can in principle be broken u sing 
magnification effect (see iBroadhurst et al.lll995ll2005|) . In 
order to make full use of this method, the unlcnscd source 
counts at a given magnitude threshold must be known 



2 



M. Bradac et al.: The combined strong and weak lensing cluster mass reconstruction method 



accurately. Given the photometric calibration uncertain- 
ties, which at the faint magnitudes one is usually deal- 
ing with easily amount to 0.1 mag, thus an uncertainty 
of ~ 10 % in the unlensed sou rce counts is typical. As 
shown bv lSchneider et all l)2000|h this level of uncertainty 
removes a great deal of the power of the magnification 
method to break the mass-sheet degeneracy. 

Several attempts have been made recently to measure 
the clust er mass profiles with weak lensi ng. However, as 
shown in IClowe fc Schneider! 1 200ll |2002() , it is extremely 
difficult to distinguish e.g. isothermal from NFW profiles 
at high significance using weak lensing data alone. The au- 
thors also conclude that these difficulties mostly arise as a 
consequence of the mass-sheet degeneracy transformation. 
Therefore additional information needs to be included, 
such as co mbining the weak lensing data with st rong lens- 
ing (see e.g. lKneib et a l. 2003, Smith ct al. 2004). Another 
example was given bv lSand et all J2004I) using combined 
strong lensing and stellar kinematics data of the domi- 
nating central galaxy. This approach offers valuable extra 
constra ints, however a detailed strong-lens modelling is re- 
quire d (jBartelmann fc MeneghettilEool IDalal fc KeetorJ 
120031) . 

In this paper we use a combined strong and weak lens- 
ing mass reconstruction to determine the mass and the 
mass distribution of clusters. We reconstruct the gravi- 
tational potential ip, since it locally determines both the 
lensing distortion (for weak lensing) as well as the deflec- 
tion (for strong lensing). Th e me thod extends the idea 
from lBartelmann et al! jl996t) and lSeitz etHl (|l998h . Its 
novel feature is that we directly include strong lensing 
information. Further, weak lensing reconstruction is ex- 
tended to the critical parts of the cluster and we include 
individual redshift information of background sources as 
well as of the source(s) being multiply imaged. This al- 
lows us to break the mass-sheet degeneracy and accurately 
measure the cluster mass and mass distribut ion. The 
method is tested using sim ulations, and in iBradac et alJ 
2004af) thereafter IPaper Ilh we apply it to the cluster RX 



shift 



J1347.5— 1145. In this paper we first briefly present the ba- 
sics of gravitational lensing in Sect. El In Sect. we give 
an outline of the reconstruction method (detailed calcu- 
lations are given in the Appendix). We test the method 
using N-body simulations, and we present the results in 
Sect. 0] The conclusions and summary are the subject of 
Sect. El 



2. Gravitational lensing preliminaries 

Throughout this paper we follow the notation of 
iBartehnann fc Schneider! l)200lh . who give a detailed ac- 
count of the gravitational lensing theory presented here. 

We start by considering a lens having a projected sur- 
face mass density £(#), where 9 denotes the (angular) 
position in the lens plane. We define the dimensionless 
surface density k(0) for a fiducial source located at a red- 



k{0) 



oo and a lens (deflector) at z = z d 
£(0) . ... c 2 D, 



where £ cr = 



4nGD d D d , 



(1) 



D oo, D dl and D d ^ are the angular diameter distances be- 
tween the observer and the source at z — > oo, the observer 
and the lens, and the lens and the source, respectively. 
The choice of scaling with E cr is motivated by the fact 
that lenses with £ > E cr (i.e. k > 1) are strong enough 
to form multiple images. In this paper, however, strong 
lensing refers to multiple imaging only, while weak lensing 
means that the lensing effect is treated as statistical in 
nature, although it is also applied to the lens region with 
n > 1, traditionally called the strong lensing regime. 
We define the deflection potential ij>(6) 

if>(8) = - I d 2 9' k(0') In \9 - 9'\ . (2) 

which is related to k via the Poisson equation 

V 2 V>(6>) = 2k(9) (3) 

Also the shear 7 = 71 + i72 and the deflection angle a are 
related to ip, where 



7i 



1p, 22) , 72 = "0,. 



(4) 



The relations and J2J| are written for the source at 
redshift z — > 00; we note, however, that they hold for 
any redshift. Since we will work with sources at different 
redshifts (both for strong as well as for weak lensing), we 
factorise the redshift dependence of the lens convergence 
k, the shear 7, and the deflection angle a by 



k{0,z) 
a(0,z) 



1 {9,z) = Z{z) 1 {9) , 



Z{z)k{0) , 

Z{z)ol{9) . (5) 
Z(z) is the so-called "cosmological weight" function: 
DooD d 



Z(z) 



Ai,ooA 



H(z - z d ) 



(6) 



where D ds and D s , are the angular diameter distances 
between the lens and the source, and the observer and 
the source at a redshift z, respectively. H(z — z d ) is the 
Heaviside step function and accounts for the fact that 
galaxies located at z < z d are not lensed. 1 

In the case of weak lensing, the information on the lens 
potential is contained in the transformation between the 
source ellipticity and image elliptic ity e. It is given as 
a fun ction of reduced shear g(9, z) fsee lSeitz fc Schneider! 
119971) 





f e-g(o,z) 




l-g*(0,z)e 


e (s) = < 






l-g(0,z)e* 




[ e* - g*{9, z) 



for \g{0,z)\ < 1 



for \g{0, z)\ > 1 



(7) 



1 To evaluate the angular diameter distances we assume the 
ACDM cosmology with Q m = 0.3, = 0.7, and Hubble con- 



stant H = 70km s" 1 Mpc" 
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where 



g{o,z) = 



(8) 



1 - Z{z)k{9) ' 

Galaxies are intrinsically elliptical, and therefore one 
cannot disentangle the effect of lensing from the intrinsic 
properties in Q using a single galaxy image. Hence, the 
weak lensing effect needs to be treated in statistical sense. 
More precisely, we can Taylor expand the expression (0 
(e.g. for the case of \g(6, z)\ < 1) and recall that | < 1 
by definition to obtain 



e{z) = 



As) 



9(0, z) 



1 + .g*(0,z)e( s ) 
(e( s >+ 5 (0,z 



E 

k=Q 



(-l) k (g*(0,z)) k (e^) k . (9) 



A similar expansion can be obtained for the case 
\g(6, z) | > 1. If we assume that the intrinsic ellipticity 

distribution has moments (e^ ) = for each k except 
k = we get the known expression for the expectation 
value of the image ellipticity at redshift z 




if \g(0,z)\ < 1 



otherwise 



(10) 



g*(0,z) 

This relation is particularly simple due to the convenient 
definition of e. In the approximation k <C 1, |7| <C 1 (thus 
\g\ <C 1) the expectation value is given by (e(z)) = 7(0, z). 

3. The cluster mass reconstruction method 

The idea of combining strong and weak lensing con- 
straints is not new, it has been prev i ously discussed by 
lAbdelsalam et al.ll|l998|) . lKneib et aDl|2003() . ISmith et alJ 
l)2004l) . and others. The method presented here, how- 
ever, has some import ant differences. For example in 



lAbdelsalam et all l)1998l) the authors reconstruct the pix- 
elized version of the surface mass density k. A similar 
method for strong le nsing constraints o nly has recently 
also been presented bv lDiego et alJl)2004fK We argue, how- 
ever, that using the potential tp is favourable, since k, 7, 
and a locally depend upon the potential ip - c.f. 0, l@J - 
and all can be quantified from the latter. 7 and a, on the 
other hand, are non-local quantities of k. In other words, 
the mass density on a finite field does not describe the 
shear and the deflection angle in this field. If a finite field 
is used, one usually employs Fourier analysis; in this case, 
7 in fact corresponds to original k plus all its periodic 
continuations. 

Further, even though not easy to implement, we de- 
cided to keep the parametrisati on of the mas s- distrib ution 
as general as po ssible. In iKneib et alJ I^OO^ and 
ISmith et alJ l|2004l) . on the other hand, the strong and 
weak lensing constraints were compared in a Bayesian ap- 
proach in the form of simple, parametrised models. In 
addition, the weak lensing constraints were not used to 
the very centre of the cluster and rcdshifts of individual 
sources were not included. 



3.1. The outline of the method 

The main idea behind the method is to describe the cluster 
mass-distribution by a fully general lens, using the values 
of the deflection potential ip on a regular grid. We then 
define a penalty function \ 2 and minimise it with respect 
to the values of iji^. The convergence K, the shear 7, and 
the deflection angle a at an arbitrary position in the field 
are obtained by finite differencing and bilinear interpola- 
tion methods. The number of grid points we use for ip/. is 
(N x + 2) x (N y + 2); the extension by one row and one 
column at each side is needed to perform the finite differ- 
encing at each inner N x x N y grid point. 
We define the x 2 -function as follows 



(ii) 



Xei'tpk) contains information from statistical weak lensing, 
whereas in Xm (V'fc ) we include the multiple imaging prop- 
erties. R(ipk) is a regularisation term multiplied by the reg- 
ularisation parameter r\. The regularisation is a function 
of the potential and disfavours small-scale fluctuations in 
the surface mass density. 

In order to find the minimum \ 2 solution, we solve 







(12) 



This is in general a non-linear set of equations, which we 
solve in an iterative manner. We linearise this system and 
in the first step we start from some trial solution, to calcu- 
late its non- linear terms (see the Appendix for details). We 
solve the corresponding equations and repeat this proce- 
dure until a convergence is achieved. Inverting the result- 
ing matrix of ~ (N K x A'y) 2 elements for finding a solution 
of the linear system is difficult in general even for grids 
with a small number of cells. However, as it turns out, 
the resulting matrix is sparse and the system is therefore 
computationally inexpensive to solve. 

The reconstruction is performed in a two-level iteration 
process, outlined in Fig. ^ We will refer to the iteration 
process mentioned above for solving the linear system of 
equations as inner-level, where steps ri\ are repeated until 
convergence of k. The outer-level iteration is performed 
for the purpose of regularisation (as described in detail 
in Sect. In order to penalise small-scale fluctuations 
in the surface mass density, we start the reconstruction 
with a coarse grid (large cell size). Then for each rii step 
we increase the number of grid points in the field and 
compare the new reconstructed k'" 2 -* with the one from 
the previous iteration k^" 2-1 ) (or with the initial input 
value for 712 = 0). The second-level iterations are per- 
formed until the final grid size is reached and convergence 
is achieved. 



3.2. Technical aspects 

In this section we will briefly describe some technical as- 
pects of how we calculate the lensing quantities k, 7, and 
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n 2 = 



! 

— | Initial coefficients of the matrix with ct(" 2 ) n( n -^ 7 (" 2 ) 



Recalculate coefficients of the matrix 



outer-level 
iteration 

repeat 112 until fi- 
nal grid size 







inner-level 




Solve the linear 


iteration 




system 


repeat n± until 






convergence 



Using ifr, calculate a ( - ni+1 '> 7 ("i+ 1 ) [ 



Increase the num- 
ber of grid points 



[interpolate qt" 2+1 ) ^("a+i) 7 (na+i) on a new g^jj 

Fig. 1. The outline of the two- level iteration process. 



a at an arbitrary position within the field from the po- 
tential ipk on a grid. 

We use the finite differencing method with 9 grid 
points to calculate n, 5 point s for 7, and 4 points for a (see 
lAbramowitz fc Stegunlll97^ . The coefficients used for n 
and 7 are given in Fig. the case of a is discussed in 
Appendix IA. 21 To evaluate k(9), 7(0), and a(0) at a po- 
sition 9 within the field, bilinear interpolation is used. 

Note, that the dimensionality of the problem is not 
(N x + 2) x (N y + 2). Because the transformation ip(9) — » 
ip(9) + ipQ + a ■ 9 leaves n and 7 invariant, the poten- 
tial needs to be fixed at th ree points (see ISeitz et all 
ll998llBartehnann et alll99fij) . These thus fix the constant 
and linear term in the invariance transformation. If this 
is not the case, a minimum in x 2 would correspond to 
a three-dimensional subspace of possible solutions. The 
choice of the three points, and the corresponding values 
of the potential are arbitrary. Although the transforma- 
tion ip(9) — > rp(8) + a -9 changes the deflection angle ex, it 
only causes a translation of the source plane, which is not 
an observable. Therefore, even in the presence of strong 
lensing, three points of the potential need to be held fixed. 

The mass-sheet degeneracy transformation of the po- 
tential is given by ip — > ip' = (1 — X)9 2 /2 + Aip. However 
since we a i m at lifting this degeneracy, in contrast to 
ISeitz et all (^998) the potential ipk, is not held fixed at an 
additional, fourth point. The dimensionality of the prob- 



lem is thus Nri 



(N x + 2)(N y 



3.3. The ^-function 

In this section we will describe contributions to the x 2 - 
function, starting with the statistical weak lensing. 



For N g galaxies with measured ellipticities e$ we define 
the x 2 as 



X, 



where 



(13) 



(14) 



Note that (e) refers to the expectation value of the ellip- 
ticity over redshift space, not to an ensemble average over 
different galaxies and is derived from (|10|l . 

In lBradac et alJ l|2004p|) 2 we argue that x 2 can g ive bi- 
ased results for lenses for which many galaxies have \g\ — 
1. It would be better to work with a log-likelihood function 
with a probability distribution that properly describes 
the distribution of observed ellipticities. Unfortunately, 
such an approach is not viable here (as will become obvi- 
ous later on). However, in general clusters do not have 
a l arge fraction of galax ies with \g\ ~ 1 and we show 
in iBradac et alJ l)2004bj) that for these lenses the x 2 ~ 
minimisation is sufficient. 

One of the major strengths of this statistical weak 
lensing reconstruction technique is the possibility to si- 
multaneously include constraints from multiple image sys- 
tems to the weak lensing data in a relatively straight- 
forward manner. The simplest approach to strong lens- 
ing is to perform the so-called "source plane" modelling; 
i.e. to minimise the projected source position difference. 
Consider a multiple image system with the source at red- 
shift z s and with A^m images located at 9 m . We define 
b m = 9 rn — Z(z s )a(9 rn ) — f3 s and the corresponding x 2 - 



We would like to point out here that there is a typo in that 
paper. All Wi = l/cr, in the text and plot labels need to be 
replaced by Wi = 1/of , which has been used in the calculations. 
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Fig. 2. The finite differencing coefficients for k (left), 71 (middle) and 72 (right). E.g. for k we use a formula including 9 
points, the multiplicative factor is given at the bottom, the individual coefficients in the circle. Thus for the middle point 
(0, 0) we get «(0, 0) = ±^ [2 ty>(-l, 1) + V(l, 1) + - 1 ) + -1)] - [1>(Q, 1) + V>(-1, 0) + V(l, 0) + ^(0, -f )] - 

#(0,0)]. 



function is given by 

xm = E b l s ~ lb ™ . ( 15 ) 

m— 1 

where S is the covariance matrix. (3 S is the average source 
position and for simplicity we calculate it using the deflec- 
tion angle information from the previous iteration ni — 1. If 
the measurement errors on image positions are distributed 
isotropically, S reduces to a diagonal matrix given by 
S = diag((Tg 1 , <7g 2 ) with = a^ 2 - °s,i and er s ,2 are 
the errors on image positions, projected onto the source 
plane. For simplicity, however, we do not perform a pro- 
jection of the error ellipse from the image plane onto the 
source plane. Instead, we keep S constant throughout the 
reconstruction. Therefore we avoided the numerical prob- 
lem of a diverging function if one of the images lies 
at the position of the critical curve for the corresponding 
redshift. 

We are aware of the fact that the a pproach we use 
is not optimal (see e.g. iKochanekl 12004^ . If only multi- 
ple imaging is used, the resulting best-fit model is bi- 
ased towards high magnification factors, since errors on 
the source plane are magnified when projected back to 
the image plane (this information we do not use). In our 
case, however, the model also needs to take into account 
the constraints from statistical (weak) lensing and there- 
fore the high magnification models are in fact discarded. 
In addition, if e.g. one considers the errors matrix in the 
image plane to be diagonal, the corresponding matrix 
in the source plane would have large off-diagonal terms. 
Throughout this paper we therefore consider the errors in 
the source plane to be isotropic, since this may in fact be 
a better approximation, as sources are on average more 
circular than their lensed images. In practice the location 
of the multiple images are usually known very accurately, 
leading to a very narrow minimum of Xm m f ue parameter 
space. In practice, multiple image constraints are satisfied 



nearly perfectly and exact values of errors on image posi- 
tions are of lesser importance. 

Since the minimisation of Xe can lead to a potential 
that reconstructs the noise in the data, the solution needs 
to be regularised. Even without measurement errors, the 
intrinsic ellipticities would still produce pronounced small- 
scale noise peaks in the final reconstruction. In addition, 
the method presented here has an intrinsic invariance if no 
multiple imaging information is used and the weak lensing 
approximation g ~ 7 applies. Namely, we can alternately 
add/subtract a constant a along diagonals of the potential 
(chess-board like structure, as sketched in Fig. |3J). This 
transformation would on the one hand not affect 7, but on 
the other it would cause a similar change (with a constant 
2a/3) in k - compare with Fig. Thus in the \g\ -C 1 
regime, where (e) = g ~ 7 these stripes would show up in 
the resulting re-map. This problem can, however, be very 
efficiently cured with regularisation. 




Fig. 3. The intrinsic invariance of the method. If we al- 
ternately add/subtract a constant a along the diagonals 
the shear 7 does not change (cf. Fig.|2J), but n changes in 
the similar way with a constant now being 2a/3 . 

Since we want to measure the cluster mass, the reg- 
ularisation should not influence breaking the mass-sheet 
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degener acy. For ex a mple, one of the possibilities consid- 
ered bv lSeitz et al. I lll998h for regularisation function was 

R = J2fj^i r \^ K \ 2 - However, as the authors mentioned, 
such regularisation would tend to flatten the profile and 
therefore affect the mass-sheet degeneracy breaking. Their 
maximum-entropy regularisation with moving prior (i.e. 
the prior in the regularisation is not kept constant, but 
adapted in the process of minimisation) does not flatten 
the profile, however it is very cumbersome to express its 
derivative in linear terms of ipk- Motivated by the success 
of moving prior in maximum-entropy regularisation, we 
choose a very simple prescription for the regularisation 
function. We start off by a relatively coarse grid, since 
if the number of grid points A^j; m is much smaller than 
the number of galaxies, the resulting reconstruction is not 
able to follow the noise pattern. In each second-level iter- 
ation step we gradually increase the number of grid points 
and compare the resulting k map k*- 712 -* with that from the 
previous iteration k^™ 2-1 ) linearly interpolated on a finer 
grid, thus 

R = E (4 2) -4 2 ~ 1] ) ■ ( 16 ) 

For the case of n% = we use an initial guess for k which 
can in practice be obtained from strong lensing, direct 
finite-field reconstruction, parametrised model fitting to 
weak lensing data, or simply set to «y = 0. This form 
of regularisation is relatively easy to implement and is in 
addition very efficient in removing the stripes (mentioned 
above) in the final reconstruction. If enough n 2 iteration 
steps are used this form of regularisation also does not 
substantially affect the ability to break the mass-sheet de- 
generacy, since the information of initial k}- q > is lost and 
different initial conditions do not bias the results (if the 
signal from lensing is high enough). 

Finally a word on the regularisation constant 
This parameter should, in theory, be set such to ensure 
X 2 /-/Vdof ~ 1. In practice, however, it is difficult to de- 
termine its optimal v alue (in the cr i tical le nsing regime). 
As outlined in iGeieer fc Schneider! l|l998|) the probabil- 
ity distribution of measured ellipticities is not a Gaussian 
and therefore the minimum value of \ 2 has no particu- 
lar meaning. In practice, setting 77 such that the resulting 
X 2 /Ndot ~ 1 (where N^oi is in our case the number of 
galaxies used for strong and weak lensing) is valid is a 
good guess for this parameter. In addition, one adjusts 
rj low enough for the method to have enough freedom to 
adapt to the information in the data and large enough for 
not allowing the solutions that follow the noise pattern. 
As a rule of thumb it is usually better to set rj high and 
increase the number of iterations and hence allowing n to 
change only slowly. Since the reconstruction is done in a 
two-level iteration and in addition multiple-image infor- 
mation is included, the method can successfully adapt to 
the data and the results are not very sensitive to the pre- 
cise value of 77. The resulting smoothness level of the mass 
maps should reflect the quality of data. The "smoothing 



scale" depends upon the combination of the grid size and 
regularisation. The final potential map should be void of 
any structures on scales smaller than the mean separation 
between galaxies used for weak lensing. We will shortly 
return to this point in Sect. 14.31 

3.4. Initial conditions 

For the purpose of the regularisation, given by we 
need to employ the initial conditions. In addition, if a re- 
alistic model is used for the initial conditions (and as trial 
solution), it is very helpful to break the internal degen- 
eracy, i.e. to distinguish galaxies that have \g\ < 1 from 
those with \g\ > 1 in the first step, and thus allow for a 
faster convergence. Breaking this degeneracy is desirable, 
since otherwise the method has difficulties in "climbing" 
over the \g\ = 1 region (especially if no multiple- image 
systems are included) . Since we use individual redshifts of 
background galaxies we do not have well-defined critical 
curves (i.e. positions in the source plane where \g\ — 1), 
as their position depends on the source redshift. In spite 
of this fact, the transition still poses a difficulty. 

In our case different initial conditions are employed. 
For the initial model k' ' we use three different scenarios: 
= (and = 0, 7W = 0) across the whole field 
(hereafter 10), taken from the best fit non-singular 
isothermal ellipsoid model NIE for the multiple image 
system described in Sect. 14.21 (hereafter IM) and a non- 
singular isothermal sphere model NIS with scaling and 
core radius being the same as in IM (hereafter IC). The 
same models are used also to obtain the initial coeffi- 
cients of the linear system (see Appendix) (for 10 we use 
71,2 = 0). These different initial models help us to ex- 
plore the effects of regularisation and the capability of the 
reconstruction method to adapt to the data. 

4. Cluster mass reconstruction from simulated 
data 

4.1. Mock catalogues 

We generate mock catalogues using a cluster from 
the hi gh-resolution N-body simulation by ISpringel et alJ 
l|200lj) . The cluster is taken from the S4 simulation (for 
details see the aforementioned paper) and was simulated 
in the framework of the ACDM cosmology with density 
parameters J7 m = 0.3 and Q\ = 0.7, shape parame- 
ter r = 0.21, the normalisation of the power spectrum 
as = 0.9, and Hubble constant Ho — 70 kms' 1 Mpc -1 . 
The cluster simulation consists of almost 20 million parti- 
cles, each with a mass of 4.68 x 10 7 M Q and a gravitational 
softening length of 0.7 hr 1 kpc. Due to the high mass res- 
olution, the surface mass density K-map can be obtained 
by directly projecting the particles (in our case along the 
z-axis) onto a 1024 2 grid (of a side length 4 Mpc) using 
the NGP (nearest gridpoint) assignment. 

In what follows we try to generate the weak and strong 
lensing data to resemble as close as possible the data on 



M. Bradac et al.: The combined strong and weak lensing cluster mass reconstruction method 



7 




Fig. 4. The gravitational lensing properties of a simulated cluster used for generating mock catalogues for statistical 
weak lensing and for the multiple image system, a) The surface mass density k, b) the absolute value of the reduced 
shear \g\, both for a source at z s — > oo, are given in gray-scale and contours. The stars in a) denote the image positions 
of a four-image system at z s — 1.76 which we use for the reconstruction. 



the cluster RX J1347.5-1145 we will use in lPaper III The 
surface mass density of the cluster is therefore scaled to 
have a sizeable region where multiple imaging is possible 
within a 3.8 x 3.8 arcmin 2 field, for sources at redshifts 
z > 1 and a cluster at z& = 0.4. The Einstein radius for 
a fiducial source at z — > oo is roughly 6>e ~ 1', giving a 
line-of-sight integrated mass within this radius of ~ 5 x 
10 14 Mq. The cut-out of the resulting k map (for z — > oo, 
thus Z{z) — 1) we use can be seen in Fig. QJ,. 

The len sing properties are calculated as described 
in detail in IBradac et alJ l)2004cj) . The Poisson equation 
for the lens potential if) - c.f. Eq. © - is solved on 
the grid in Fourier space with a DFT (Discrete Fourier 
Tr ansformation) method us ing the FFTW library written 
bv iFrigo fc Johnson! l|l998() . The two components of the 
shear 71.2 and the deflection angle a are obtained by fi- 
nite differencing methods applied to the potential tp. These 
data are then used to generate the weak lensing catalogues 
as well as the multiple image systems. The absolute value 
of the reduced shear (again for a source with z — » 00) is 
shown in Fig. 0Jd. 

The weak lensing data are obtained by placing N g 
galaxies on a 3.8 x 3.8 arcmin 2 field. We have simulated 
two different catalogues, one with N g = 148 galaxies with 
positions corresponding to those from R-band weak lens- 
ing data of the cluster RX J1347.5— 1145 and one with 
N„ = 210 g alaxies corresponding to the I-band data used 
in lPaper il In this way we simulate the effects of "holes", 
resulting from cluster obscuration and bright stars in the 
field. 

The intrinsic ellipticities e s are drawn from a Gaussian 
distribution, each component is characterised by a — 
o> = 0.2. We use the same redshifts as those measured in 
the R and I-band data, respectively. The catalogues have 



average redshifts for background sources of (z\) = 1.18 
and (zr) = 1.14. The corresponding cosmological weights 
are evaluated assuming the ACDM cosmology (the same 
parameters are used as for the cluster simulations). 

The measurement errors e crr on the observed elliptic- 
ities are drawn from a Gaussian distribution with a = 
Ccrr = 0.1 (each component) and added to the lensed el- 
lipticities. We considered also the measurement errors on 
the redshifts of the galaxies to simulate the use of pho- 
tometric redshifts. Th ese have tj zcrr = 0.1 (1 + z,) (see 
iBolzonella et al J I2OO0I) : in adding the errors we ensured 
that the resulting redshifts are always positive. We have 
also simulated the presence of outliers in the redshift 
distribution, 10% of our background sources (chosen at 
random) are considered outliers, for these we randomly 
add/subtract Az — 1 to their redshifts (which already in- 
clude random errors). The lensed ellipticities are obtained 
using @ and interpolating the quantities k, and 7 at the 
galaxy position using bilinear interpolation, considering 
the redshifts including errors. In contrast, for the purpose 
of reconstruction we then consider galaxies to be at their 
"original" redshifts (thus equal to the observed redshifts 
in the data of RX J1347.5-1145). 

4.2. Multiple imaging 

To obtain a four-image system fro m the simulation we us e 
the method described in det ail in IBradac et alJ l|2004ch . 
With the MNEWT routine from lPress et al we solve 

the lens equation for a given source position inside the 
asteroid caustic. The source is assumed at a redshift of 
z s = 1.76. Once we have the image positions, their magni- 
fications are calculated and the fifth image is eliminated, 
since it is usually too dim and would not be observable. 
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The errors on image position s can be c onservatively es- 
timated (for the data we use in IPaper it to ~ 0'.'3. Since 
we need errors in the source plane, we set them by a fac- 
tor of five smaller (in agreement with the average magni- 
fication factor for this system), i.e. er sm = 0'.'06 for both 
coordinates (see discussion in Sect. I3.3L 

We also use this system to obtain one of the mod- 
els, needed for the purpose of strong and weak lensing re- 
construction. We perform image plane minimisation and 
fit an NIE model ijKormann et al.lll994|) . given by 

K (e') = b ° , (l?) 

where 0' is calculated w.r.t. the semi-major axis of the 
cluster surface mass density. We allow the scaling bo, 
ellipticity |e g | and position angle (f> g to vary. The best 
fit model for this system has values of {bo, |e g | , <fi g } = 
{0'.97, 0.30, 1.01}. We fix the co re radius r r to (XI. For 
model fitting we use C-minuit l|james fc Rooslll975l) . a 
routine which is part of the CERN Program Library. 

4.3. Weak lensing mass-reconstruction using simulated 
data 

The mock catalogues are used to test the performance 
of the reconstruction method. We obtain the solution of 
the linear system with the UMFPAC K routine for solvin g 
asymmetric sparse linear systems l|Davis fc Dufj 119 99). 
We perform 30 second-level iterations, each time increas- 
ing 7V X and N y by one, starting with an initial 20 x 20 
grid (since we use regularisation from the first iteration, 
we do not need to start the reconstruction with a very 
course grid, where the resulting matrix is not sparse and 
a different routine should be applied to find a good so- 
lution). The results of the reconstructions are shown in 
Fig.0for N g — 210 galaxies distributed in the same man- 
ner as the I-band data and N g = 148 galaxies d istributed 
as the R-band data of the cluster described in IPaperTB . 
The initial regularisation parameter was set to rj — 400 
for I-band rj — 200 for R-band data and adapted in each 
step to ensure xV-^dof ~ 1- Two different initial values 
for the regularisation parameter are used since the num- 
bers of galaxies in the two bands are different. It is very 
comforting to observe that the reconstructed maps do not 
depend crucially on the initial K-model we use. 

From the reconstructed maps we estimate the mass 
within the 1'.5 radius from the centre of the cluster (for 
a redshift z& — 0.4 this corresponds to 340/i _1 kpc) pro- 
jected along the line-of-sight. For this purpose we generate 
10 mock catalogues for each band and did the reconstruc- 
tion again with the three different initial models. We list 
the resulting average mass obtained from the catalogues 
in Table ^ for both the I- and R-band mock catalogues. 
All the mass estimates are similar; note, however, that the 
galaxy catalogues following the I- and R-band data have 
galaxies partly in common and the errors are therefore cor- 
related. Wc find the enclosed mass of the simulated cluster 



to be (1.0±0.1) x 1O 15 M , which is very close to the input 
value of M s {< 340/1" 1 kpc) = 0.99 x 10 15 A/ o . The l-o er- 
ror is estimated from the variance of mass determinations 
from different mock catalogues. 

The results show that our method is effectively able 
to break the mass-sheet degeneracy and is, as a conse- 
quence, very efficient in reproducing the cluster mass also 
at radii significantly larger than the Einstein radius of the 
cluster. It is also very encouraging that the results are 
nearly independent of the initial k' ' used for the regu- 
larisation. Note that a single multiple-image system does 
not by itself break this degeneracy, we would need at least 
two different redshift multiple image systems to break the 
mass-sheet degeneracy with strong lensing data alone. In 
such a case the strong lensing gives constraints on the mass 
enclosed within the Einstein radius for a given source red- 
shift and since the critical curves depend on the source 
redshift, we can constrain mass at two different radii and 
the degeneracy is broken. The combiniation of weak and 
strong lensing is the more powerful, the more different the 
redshift of the source of the multiple images is from the 
median redshift of the galaxies from which the weak lens- 
ing measurements are obtained, and the less symmetric 
the arrangement of multiple images is w.r.t. the center of 
the cluster. 

Unfortunately we can not resolve both clumps present 
in the simulations. This is due to the fact that the num- 
ber density of background sources is low and the inter- 
nal smoothing scale (i.e. the average distance between 
two source galaxies) is large; with a number density of 
~ 100 arcmin -2 the clumps can be easily resolved. 

We have also performed additional reconstructions in 
which we multiplied the original values of k of the simu- 
lated cluster by 0.75 and 1.25. This enables us to confirm 
that the agreement between input mass and reconstructed 
mass is not just accidental. We have generated new mul- 
tiple image systems and new mock catalogues as before. 
We do not, however, perform a new strong-lensing recon- 
struction, for k' ' we intentionally use the same (i.e. in 
this case "wrong") initial conditions as before. The old 
IM model would not fit the image positions any longer, 
since they have changed with the scaling of k. The recon- 
structed masses of the increased k simulation are in good 
agreement with the input values. The differences between 
different models are comparable (slightly smaller) to the 
ones shown in Table ^ For the lower k simulation, the 
reconstructed values are on average the same as the in- 
put value, however the scatter is larger. This is expected, 
since the lens in this case is weaker and the breaking of 
the mass-sheet degeneracy is difficult in this case (with 
the quality of data used here) . 

As an additional test we also consider a redshift distri- 
bution with (z) = 1.6 for the sources used for weak lensing 
and regenerate the mock catalogues. The accuracy of the 
determination of the enclosed mass increases. However, 
more importantly, we also better reconstruct the shape of 
the mass distribution, since high-rcdshift galaxies (when 
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their shape is measured reliably) contribute most to the 
signal and improve the accuracy of the reconstruction. 

5. Conclusions 

In this paper we dev elop a new method based on 
iBartelmann et all l|l996t) to perform a combined weak and 
strong lensing cluster mass reconstruction. The particular 
strength of this method is that we extend the weak lens- 
ing analysis to the critical parts of the cluster. In turn, 
this enables us to directly include multiple imaging in- 
formation to the reconstruction. Strong and weak lensing 
reconstructions are performed on the same scales, in con- 
trast to similar methods proposed in the past where weak 
lensing information on much larger radii than the Einstein 
radius 6 k was combined w ith strong lensing information 
(see e.g. lKneib et~al1l2003|) . 

We test the performance of the method on simulated 
data and conclude that if a quadruply imaged system com- 
bined with weak lensing data and individual photometric 
redshifts is used, the method can very successfully recon- 
struct the cluster mass distribution. With a relatively low 
number density of background galaxies, 15 arcmin -2 , we 
are effectively able to reproduce the main properties of the 
simulated cluster. In addition, with larger number densi- 
ties ~ 100 arcmin -2 of background sources, accessible by 
HST, the substructures in the cluster can be resolved and 
the mass determination further improved. 

We determine the enclosed mass within 340ft" 1 kpc of 
the simulated cluster to be (1.0 ± 0.1) x 10 15 M Q , which 
is very close to the input value of M s (< 340ft -1 kpc) = 
0.99 x 1O 15 M . We have shown, that with the data qual- 
ity we use we are effectively able to break the mass- 
sheet degeneracy and therefore obtain the mass and mass- 



Table 1. Reconstructed cluster mass within a cylinder of 
340ft. -1 kpc radius around the cluster centre from simu- 
lations of mock catalogues resembling I-band (left) and 
R-band (right) weak lensing data and one 4-image sys- 
tem. Three different initial conditions are used. We use 
the best fit model from the multiple image system IM 
(see Sect. I4.2|l . the IC model (NIS with same scaling and 
core radius as IM) and 10 with — on all grid points. 
In the last line the input mass M s from the simulation is 
given. The variance of the mass estimate is given and in 
brackets we give for comparison the velocity dispersion of 
an SIS having the same enclosed mass within 340ft. -1 kpc. 







Mi, s 
[1O 15 M ] 


[<TI,SIS] 

[kms -1 ] 


M R , B 
[1O 15 M ] 


[CTR.SIS] 

[kms" 1 ] 


IM 
IC 
10 


1.07 ±0.02 
1.04 ±0.02 
0.88 ± 0.03 


[1720] 
[1700] 
[1560] 


1.04 ± 0.02 
1.00 ±0.02 
0.82 ±0.03 


[1710] 
[1670] 
[1510] 




M s (< 340ft" 1 kpc) 


0.99 


[1670] 



distribution estimates without prior assumptions on the 
lensing potential. 

In addition, the reconstruction algorithm can be im- 
proved in many ways. First, we use for the multiply im- 
aged system only the information of the image positions. 
The reconstruction method can, however, be modified to 
include the morphological information of each extended 
source. Instead of using a regular grid, one would have to 
use adaptive grids and decrease the cell sizes around each 
of these images. This will be a subject of a future work. 
Second, the photometric redshift determination does not 
only give the most likely redshift given the magnitudes in 
different filters, but also the probability distribution for 
the redshift. This information can be included in the re- 
construction. In addition, source galaxies without redshift 
information can be included and different regularisation 
schemes can be considered. 

Finally, the slight dependence on the initial conditions 
is getting weaker the higher the number density of back- 
ground galaxies and/or multiple image systems are. In ad- 
dition, it is of advantage to have a large spread in the red- 
shift efficiency factors Z of the background galaxies. For 
example, deep ACS images of clusters with a usable num- 
ber density of n ~ 120 arcmin" 2 , or future observations 
with the James Webb Space Telescope will most likely 
make the depe ndence on the initial conditions negligible. 

In lPaperU we will show the application of this method 
on the cluster RX J1347.5— 1145 and confirm that a com- 
bination of strong and weak lensing offers a unique tool 
to pin down the masses of galaxy-clusters as well as their 
mass distributions. 
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Appendix A: The linear problem for ^ 

In this section we present details of the method outlined 
in Sect. 

We aim to solve the equation 



dipk 



v- 



= 



(A.l) 



This is in general a non-linear system of equations. We try 
to solve it in an iterative way by linearising the equation 
in terms of i/'fe and keeping the non-linear terms fixed at 
each iteration step. The resulting system is then written 
in the form 



(A.2) 
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where the matrix Bjk and vector Vj contain the contribu- A. 2. The strong-lensing term 



tions from the non-linear part. In the following sections 
we will describe the contributions to i|A.l[l in turn. 



A.l. The weak lensing analysis 

The Xe f° r the weak lensing case is given in JT3J. From 
now on we consider in detail only the \g\ < 1 case; for 
\g\ > 1, the calculations are done in the same fashion. 
First we plug into (fH5|l the expectation value of observed 
ellipticities (i.e. the reduced shear g) obtaining 



E 

»=i 



X-Zk 



E 



\e-ZeK-Z-fY 
(1-Zk) 2 ct 2 



(A.3) 

where k, 7, and a depend on 6i only and Z depends on 
the redshift of the i-th source. Note that for simplicity we 
omit the index i to e, k, 7, a and Z , although these quanti- 
ties are different for every galaxy. As described in Sect. 13.21 
using finite differencing and bilinear interpolation, we can 
write k and 7 at each galaxy position as a linear com- 
bination of ipk ■ This is expressed in the following matrix 
notation 



7i(#i) = Gil'i/ik , 72 = Gil Vfc ) K (9i) = fciki^k , 

(A.4) 

where the matrices G^ , G^ , and ICik are composed of 
numerical factors described in Sect. 13.21 Now we consider 
the denominator of (|A. 3|> it< = (1 — Zk) 2 a 2 fixed at each 
iteration step (the subscript < denotes the \g\ < 1 case) 
and differentiate the following term of ljA.31) 



1 d \e - Zen - Zrf 
Z 



+ (e 2 - Ze 2 K - Z72) e 2 



Ok 
'dtfjk 



^72 

dipk 



V'fe 



+ e< J + (e? + 4) *< 



,(2) 



(A.5) 



where ei and e 2 are the two components of the measured 
ellipticity of galaxy i (again omitting the index and we 
divided Eq. IjA.lH by two for simplicity). We sum over all 
galaxies used for the weak lensing analysis and obtain a 
linear problem for ipk at each iteration step. The same 
approach can be used for the \g\ > 1 case, where cr2 is 



kept constant and is given by ai 



Following the prescription from the previous section we 
now write the deflection angle in a matrix form 



ai(0 m ) = X>£V* , a 2 (0 m ) = Vf^ k 



(A.6) 



Both matrices give the finite differencing form for the gra- 
dient of the potential, in particular we use the central dif- 
ferencing formula, i.e. ai(Q,0) = ^(^(1,0) — 1,0)) 
and a a (0, 0) = 33 (^(0, 1) - V(0, -1)). 

The \ 2 contribution to strong lensing is given in (|15JI . 
The source position f3 s is kept constant at every iteration 
step, and is evaluated using the deflection angle informa- 
tion o^ 711-1 ) from the previous iteration 



, N M 
rn—i 



(A.l) 



We differentiate the following term in %m (f° r tne x i~ 
coordinate) 

1 d {9 m , x -Z{z a )ax{e m )-^,xf 



s,m(l) 



m ,x - - Z(z a )VY'V\^ k 



(A.8 



s,m(l) 



The expression for the a; 2 -coordinate is obtained by ex- 
changing 1 — > 2. After summation of both terms over all 
images m we get a set of equations which are linear in ip k 
and can be readily included in 1A.2|I . In principle we could 
also use a(tpk) in (|A.7J| . the expression dxli/dipk would 
remain linear in tpk- However since the first approach is 
computationally much simpler, we use the former. 

A.3. The final result 

In the previous section we describe how we linearise the 
contributions of weak and strong lensing, now we can write 
the coefficients in the equation l|A.2j) . Note that the con- 
tribution of the regularisation term (with ^-contribution 
given in ]l<j\i ) is already linear in ?pk and therefore the full 
matrix Bjk is given in the form 



E 

i=l 



(2)^(2) 
ik 



+a 13 (i) [GlflQk + G^lQ. 

+a33(«) {fcijfcik)] 



N M 



(2)r,(2) 

mk 



m— 1 



(A.9) 



where the sums over i, g and m denote summation over 
all galaxies with ellipticity measurements, all grid points, 
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Fig. 5. K-maps obtained from combined strong and weak lensing reconstruction of the simulated data. Left panels 
show the reconstructions using N g = 210 galaxies distributed in the same manner as the I-band, while for the right 
pan els we use N g = 148 galaxies distributed in the same manner as the R-band weak lensing data for RX J1347.5— 1145 
fsee lPaper iJ) . The galaxies have been lensed by an N-body simulated cluster. Different initial conditions are used for 
the reconstruction. In al-a2) we use best fit model from the multiple image system IM (see Sect. I4.2|l in bl-b2) we 
use the IC model, an NIS model with the same scaling and core radius as IM and in cl-c2) we use 10, i.e. = on 
all grid points. The positions of the cluster centre and two major subclumps are plotted as white circles. 



